systemPipeR: Workflow Design and Reporting Environment
Last update: 23 April, 2022
Package
systemPipeR 2.1.27
Source code download: [ .R ]
0.0.1 Introduction
systemPipeR
provides flexible utilities for designing, building and running automated end-to-end
analysis workflows for a wide range of research applications, including
next-generation sequencing (NGS) experiments, such as RNA-Seq, ChIP-Seq,
VAR-Seq and Ribo-Seq (H Backman and Girke 2016). Important features include a uniform
workflow interface across different data analysis applications, automated
report generation, and support for running both R and command-line software,
such as NGS aligners or peak/variant callers, on local computers or compute
clusters (Figure 1). The latter supports interactive job submissions and batch
submissions to queuing systems of clusters.
systemPipeR has been designed to improve the reproducibility of large-scale data analysis
projects while substantially reducing the time it takes to analyze complex omics
data sets. Its unique features allow the creation of a uniform workflow interface
and management system that allows the users to run selected steps, customize,
and design entirely new workflows. Also, the package features take advantage of
central community S4 classes of the Bioconductor ecosystem and command-line-based
software support.
The main motivation and advantages of using systemPipeR for complex data analysis tasks are:
- Facilitates the design of complex workflows involving multiple R/Bioconductor packages
- Common workflow interface for different applications
- Makes analysis with Bioconductor utilities more accessible to new users
- Simplifies usage of command-line software from within R
- Reduces the complexity of using compute clusters for R and command-line software
- Accelerates runtime of workflows via parallelization on computer systems with multiple CPU cores and/or multiple compute nodes
- Improves reproducibility by automating analyses and generation of analysis reports
Figure 1: Relevant features in systemPipeR.
Workflow design concepts are illustrated under (A). Examples of systemPipeR's visualization functionalities are given under (B).
A central concept for designing workflows within the systemPipeR environment
is the use of workflow management containers. Workflow management containers allow
the automation of design, build, run and scale different steps and tools in data analysis.
systemPipeR adopted the widely used community standard
Common Workflow Language (CWL) (Amstutz et al. 2016)
for describing parameters analysis workflows in a generic and reproducible manner,
introducing SYSargsList workflow control class (see Figure 2).
Using this community standard in systemPipeR has many advantages. For instance,
the integration of CWL allows running systemPipeR workflows from a single
specification instance either entirely from within R, from various command-line
wrappers (e.g., cwl-runner) or from other languages (, e.g., Bash or Python).
systemPipeR includes support for both command-line and R/Bioconductor software
as well as resources for containerization, parallel evaluations on computer clusters
along with the automated generation of interactive analysis reports.
Figure 2: Overview of systemPipeR workflows management instances. A) A
typical analysis workflow requires multiple single software (red), as well the
description of the input data (green), and the expected outfiles and reports
analysis (purple). B) systemPipeR provides multiple utilities to design and
build a workflow, allowing multi-instance, integration of R code and command-line
software, a simple and efficient annotation system, that allows automatic control
of the input and output data, and multiple resources to manage the entire workflow.
c) The execution of the analysis is independent of the design and build,
enabling portability and more meaningful code sharing. systemPipeR provides
options to execute a single step or multi-steps in a compute session environment,
enabling scalability, resources to re-run one or multi-steps, checkpoints, and
generation of execution reports that can track parameters and versions, providing
transparency and data provenance.
An important feature of systemPipeR's CWL interface is that it provides two
options to run command-line tools and workflows based on CWL. First, one can
run CWL in its native way via an R-based wrapper utility for cwl-runner or
cwl-tools (CWL-based approach). Second, one can run workflows using CWL’s
command-line and workflow instructions from within R (R-based approach). In the
latter case the same CWL workflow definition files (e.g. *.cwl and *.yml)
are used but rendered and executed entirely with R functions defined by
systemPipeR, and thus use CWL mainly as a command-line and workflow
definition format rather than execution software to run workflows. In this regard
systemPipeR also provides several convenience functions that are useful for
designing and debugging workflows, such as a command-line rendering function to
retrieve the exact command-line strings for each data set and processing step
prior to running a command-line.
This overview introduces the design of a new CWL S4 class in systemPipeR,
as well as the custom command-line interface, combined with the overview of all
the common analysis steps of NGS experiments.
0.0.2 Workflow Management with SYSargsList
systemPipeR allows creation (multi-step analyses) and execution of workflow
entirely for R, with control, flexibility, and scalability of all processes.
Furthermore, the workflow execution can be integrated with compute clusters
from R, accelerating results acquisition.
The flexibility of systemPipeR's new interface workflow management class is
the driving factor behind the use of as many steps necessary for the analysis
as well as the connection between command-line- or R-based software. The
connectivity among all workflow steps is achieved by the SYSargsList workflow
management class.
SYSargsList S4 class is a list-like container where each instance stores all the
input/output paths and parameter components required for a particular data
analysis step (see Figure 3).
The SYSargsList constructor function will generate the instances, using as data
input initial targets files, as well as two-parameter files (for details, see below).
When running preconfigured workflows, the only input the user needs to provide
is the initial targets file containing the paths to the input files (e.g., FASTQ)
along with unique sample labels. Subsequent targets instances are created
automatically, based on the connectivity establish between the steps. The
parameters required for running command-line software is provided by the
parameter (*.cwl and *.yml)) files described below.
The class store one or multiple steps, allowing central control for running, checking status, and monitor complex workflows from start to finish. This design enhances the systemPipeR workflow framework with a generalized, flexible, and robust design.
Figure 3: Workflow steps with input/output file operations are controlled by
SYSargsList objects. Each SYSargs2 instance is constructed from one targets
and two param files or LineWise instance for R code-based. The only input provided by the user is the
initial targets file. Subsequent targets instances are created automatically, from the previous
output files. Any number of predefined or custom workflow steps are supported. One
or many SYSargsList objects are organized in an SYSargsList container.
0.1 Getting Started
0.1.1 Installation
The R software for running
systemPipeR
can be downloaded from CRAN. The
systemPipeR environment can be installed from the R console using the
BiocManager::install
command. The associated data package
systemPipeRdata
can be installed the same way. The latter is a helper package for generating
systemPipeR workflow environments with a single command containing all
parameter files and sample data required to quickly test and run workflows.
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("systemPipeR")
BiocManager::install("systemPipeRdata")Please note that if you desire to use a third-party command line tool, the particular tool and dependencies need to be installed and executable. See details.
0.1.2 Loading package and documentation
library("systemPipeR") # Loads the package
library(help = "systemPipeR") # Lists package info
vignette("systemPipeR") # Opens vignette0.1.3 Load sample data and workflow templates
The mini sample FASTQ files used by this overview vignette as well as the
associated workflow reporting vignettes can be loaded via the
systemPipeRdata package as shown below. The chosen data set
SRP010938 obtains 18
paired-end (PE) read sets from Arabidposis thaliana (Howard et al. 2013). To
minimize processing time during testing, each FASTQ file has been subsetted to
90,000-100,000 randomly sampled PE reads that map to the first 100,000
nucleotides of each chromosome of the A. thalina genome. The corresponding
reference genome sequence (FASTA) and its GFF annotation files (provided in the
same download) have been truncated accordingly. This way the entire test sample
data set requires less than 200MB disk storage space. A PE read set has been
chosen for this test data set for flexibility, because it can be used for
testing both types of analysis routines requiring either SE (single-end) reads
or PE reads.
The following generates a fully populated systemPipeR workflow environment
(here for RNA-Seq) in the current working directory of an R session. At this time
the package includes workflow templates for RNA-Seq, ChIP-Seq, VAR-Seq, and Ribo-Seq.
Templates for additional NGS applications will be provided in the future.
sytemPipeRdata::genWorkenvir(workflow = "rnaseq")
setwd("rnaseq")0.1.4 Project structure
The working environment of the sample data loaded in the previous step contains the following pre-configured directory structure (Figure 4). Directory names are indicated in green. Users can change this structure as needed, but need to adjust the code in their workflows accordingly.
- workflow/ (e.g. rnaseq/)
- This is the root directory of the R session running the workflow.
- Run script ( *.Rmd) and sample annotation (targets.txt) files are located here.
- Note, this directory can have any name (e.g. rnaseq, varseq). Changing its name does not require any modifications in the run script(s).
- Important subdirectories:
- param/
- param/cwl/: This subdirectory stores all the CWL parameter files. To organize workflows, each can have its own subdirectory, where all
CWL paramandinput.ymlfiles need to be in the same subdirectory.
- param/cwl/: This subdirectory stores all the CWL parameter files. To organize workflows, each can have its own subdirectory, where all
- data/
- FASTQ files
- FASTA file of reference (e.g. reference genome)
- Annotation files
- etc.
- results/
- Analysis results are usually written to this directory, including: alignment, variant and peak files (BAM, VCF, BED); tabular result files; and image/plot files.
- Note, the user has the option to organize results files for a given sample and analysis step in a separate subdirectory.
- param/
Figure 4: systemPipeR’s preconfigured directory structure.
The following parameter files are included in each workflow template:
targets.txt: initial one provided by user; downstreamtargets_*.txtfiles are generated automatically*.param/cwl: defines parameter for input/output file operations, e.g.:hisat2/hisat2-mapping-se.cwlhisat2/hisat2-mapping-se.yml
- Configuration files for computer cluster environments (skip on single machines):
.batchtools.conf.R: defines the type of scheduler forbatchtoolspointing to template file of cluster, and located in user’s home directory*.tmpl: specifies parameters of scheduler used by a system, e.g. Torque, SGE, Slurm, etc.
0.1.5 Structure of targets file
The targets file defines all input files (e.g. FASTQ, BAM, BCF) and
sample comparisons of an analysis workflow. The following shows the format of a
sample targets file included in the package. It also can be viewed and
downloaded from systemPipeR’s GitHub repository
here.
In a target file with a single type of input files, here FASTQ files of
single-end (SE) reads, the first three columns are mandatory including their
column names, while it is four mandatory columns for FASTQ files of PE reads.
All subsequent columns are optional and any number of additional columns can be
added as needed. The columns in targets files are expected to be tab separated (TSV format).
The SampleName column contains usually short labels for
referencing samples (here FASTQ files) across many workflow steps (e.g.
plots and column titles). Importantly, the labels used in the SampleName
column need to be unique, while technical or biological replicates are
indicated by duplicated values under the Factor column. For readability
and transparency, it is useful to use here a short, consistent and informative
syntax for naming samples and replicates. To avoid problems with other
packages or external software, it is recommended to use the basic naming rules
for R objects and their components as outlined here.
This is important since the values used under the SampleName and Factor
columns are intended to be used as labels for naming columns or plotting features
in downstream analysis steps.
Users should note here, the usage of targets files is optional when using systemPipeR’s new CWL interface. They can be replaced by a standard YAML input file used by CWL. Since for organizing experimental variables targets files are extremely useful and user-friendly. Thus, we encourage users to keep using them.
0.1.5.1 Structure of targets file for single-end (SE) samples
library(systemPipeR)
targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR")
showDF(read.delim(targetspath, comment.char = "#"))To work with custom data, users need to generate a targets file containing
the paths to their own FASTQ files and then provide under targetspath the
path to the corresponding targets file.
0.1.5.2 Structure of targets file for paired-end (PE) samples
For paired-end (PE) samples, the structure of the targets file is similar, where
users need to provide two FASTQ path columns: FileName1 and FileName2
with the paths to the PE FASTQ files.
targetspath <- system.file("extdata", "targetsPE.txt", package = "systemPipeR")
showDF(read.delim(targetspath, comment.char = "#"))0.1.5.3 Sample comparisons
Sample comparisons are defined in the header lines of the targets file
starting with ‘# <CMP>.’
readLines(targetspath)[1:4]## [1] "# Project ID: Arabidopsis - Pseudomonas alternative splicing study (SRA: SRP010938; PMID: 24098335)"
## [2] "# The following line(s) allow to specify the contrasts needed for comparative analyses, such as DEG identification. All possible comparisons can be specified with 'CMPset: ALL'."
## [3] "# <CMP> CMPset1: M1-A1, M1-V1, A1-V1, M6-A6, M6-V6, A6-V6, M12-A12, M12-V12, A12-V12"
## [4] "# <CMP> CMPset2: ALL"The function readComp imports the comparison information and stores it in a
list. Alternatively, readComp can obtain the comparison information from
the corresponding SYSargsList object (see below). Note, these header lines are
optional. They are mainly useful for controlling comparative analyses according
to certain biological expectations, such as identifying differentially expressed
genes in RNA-Seq experiments based on simple pair-wise comparisons.
readComp(file = targetspath, format = "vector", delim = "-")## $CMPset1
## [1] "M1-A1" "M1-V1" "A1-V1" "M6-A6" "M6-V6" "A6-V6" "M12-A12"
## [8] "M12-V12" "A12-V12"
##
## $CMPset2
## [1] "M1-A1" "M1-V1" "M1-M6" "M1-A6" "M1-V6" "M1-M12" "M1-A12"
## [8] "M1-V12" "A1-V1" "A1-M6" "A1-A6" "A1-V6" "A1-M12" "A1-A12"
## [15] "A1-V12" "V1-M6" "V1-A6" "V1-V6" "V1-M12" "V1-A12" "V1-V12"
## [22] "M6-A6" "M6-V6" "M6-M12" "M6-A12" "M6-V12" "A6-V6" "A6-M12"
## [29] "A6-A12" "A6-V12" "V6-M12" "V6-A12" "V6-V12" "M12-A12" "M12-V12"
## [36] "A12-V12"0.2 Project initialization
To create a workflow within systemPipeR, we can start by defining an empty
container and checking the directory structure:
sal <- SPRproject()## Creating directory '/home/dcassol/z_tests/GEN242_toPush/content/en/tutorials/systempiper/rnaseq/.SPRproject'
## Creating file '/home/dcassol/z_tests/GEN242_toPush/content/en/tutorials/systempiper/rnaseq/.SPRproject/SYSargsList.yml'Internally, SPRproject function will create a hidden folder called .SPRproject,
by default, to store all the log files.
In this stage, the object sal is a empty container, except for the project
information. The project information can be accessed by the projectInfo method:
sal## Instance of 'SYSargsList':
## No workflow steps addedprojectInfo(sal)## $project
## [1] "/home/dcassol/z_tests/GEN242_toPush/content/en/tutorials/systempiper/rnaseq"
##
## $data
## [1] "data"
##
## $param
## [1] "param"
##
## $results
## [1] "results"
##
## $logsDir
## [1] ".SPRproject"
##
## $sysargslist
## [1] ".SPRproject/SYSargsList.yml"Also, the length function will return how many steps this workflow contains,
and in this case, it is empty, as follow:
length(sal)## [1] 00.3 Workflow Design
systemPipeR workflows can be designed and built from start to finish with a
single command, importing from an R Markdown file or stepwise in interactive mode
from the R console.
In the next section, we will demonstrate how to build the workflow in an interactive mode, and in the following section, we will show how to build from an R Markdown file.
New workflows are constructed, or existing ones modified, by connecting each step
via appendStep method. Each SYSargsList instance contains instructions needed
for processing a set of input files with a specific command-line and the paths to
the corresponding outfiles generated.
The constructor function Linewise is used to build the R code-based step.
0.3.1 Build workflow interactive
This tutorial will demonstrate how to create a step in the workflow for running the short read aligner HISAT2 (Kim, Langmead, and Salzberg 2015).
The constructor function renders the proper command-line strings for each sample
and software tool, appending a new step in the SYSargsList object defined in the
previous step. For that, the SYSargsList function requires data from three input files:
- CWL command-line specification file (`wf_file` argument);
- Input variables (`input_file` argument);
- Targets file (`targets` argument).In CWL, files with the extension .cwl define the parameters of a chosen
command-line step or workflow, while files with the extension .yml define
the input variables of command-line steps. Note, input variables provided
by a targetsPE file can be passed on to a SYSargsList instance via the inputvars
argument of the SYSargsList function.
appendStep(sal) <- SYSargsList(step_name = "hisat2_mapping", dir = TRUE, targets = "targetsPE.txt",
wf_file = "workflow-hisat2/workflow_hisat2-pe.cwl", input_file = "workflow-hisat2/workflow_hisat2-pe.yml",
dir_path = "param/cwl", inputvars = c(FileName1 = "_FASTQ_PATH1_", FileName2 = "_FASTQ_PATH2_",
SampleName = "_SampleName_"))For a overview of the workflow, we can check the object as follows:
sal## Instance of 'SYSargsList':
## WF Steps:
## 1. hisat2_mapping --> Status: Pending
## Total Files: 72 | Existing: 0 | Missing: 72
## 1.1. hisat2
## cmdlist: 18 | Pending: 18
## 1.2. samtools-view
## cmdlist: 18 | Pending: 18
## 1.3. samtools-sort
## cmdlist: 18 | Pending: 18
## 1.4. samtools-index
## cmdlist: 18 | Pending: 18
## Note that the workflow status is Pending, which means the workflow object is rendered in R; however, we did not execute the workflow yet.
Several accessor methods are available to explore the SYSargsList object.
Of particular interest is the cmdlist() method. It constructs the system
commands for running command-line software as specified by a given .cwl
file combined with the paths to the input samples (e.g. FASTQ files) provided
by a targets file. The example below shows the cmdlist() output for
running HISAT2 on the first PE read sample. Evaluating the output of
cmdlist() can be very helpful for designing and debugging .cwl files
of new command-line software or changing the parameter settings of existing
ones.
For more details about the command-line rendered for each sample, it can be checked as follows:
cmdlist(sal, step = "hisat2_mapping", targets = 1)## $hisat2_mapping
## $hisat2_mapping$M1A
## $hisat2_mapping$M1A$hisat2
## [1] "hisat2 -S ./results/M1A.sam -x ./data/tair10.fasta -k 1 --min-intronlen 30 --max-intronlen 3000 -1 ./data/SRR446027_1.fastq.gz -2 ./data/SRR446027_2.fastq.gz --threads 4"
##
## $hisat2_mapping$M1A$`samtools-view`
## [1] "samtools view -bS -o ./results/M1A.bam ./results/M1A.sam "
##
## $hisat2_mapping$M1A$`samtools-sort`
## [1] "samtools sort -o ./results/M1A.sorted.bam ./results/M1A.bam -@ 4"
##
## $hisat2_mapping$M1A$`samtools-index`
## [1] "samtools index -b results/M1A.sorted.bam results/M1A.sorted.bam.bai ./results/M1A.sorted.bam "The outfiles components of SYSargsList define the expected output files for
each step in the workflow; some of which are the input for the next workflow step,
here next SYSargsList instance (see Figure 3).
outfiles(sal)## $hisat2_mapping
## DataFrame with 18 rows and 4 columns
## hisat2_sam samtools_bam samtools_sort_bam
## <character> <character> <character>
## M1A ./results/M1A.sam ./results/M1A.bam ./results/M1A.sorted..
## M1B ./results/M1B.sam ./results/M1B.bam ./results/M1B.sorted..
## A1A ./results/A1A.sam ./results/A1A.bam ./results/A1A.sorted..
## A1B ./results/A1B.sam ./results/A1B.bam ./results/A1B.sorted..
## V1A ./results/V1A.sam ./results/V1A.bam ./results/V1A.sorted..
## ... ... ... ...
## M12B ./results/M12B.sam ./results/M12B.bam ./results/M12B.sorte..
## A12A ./results/A12A.sam ./results/A12A.bam ./results/A12A.sorte..
## A12B ./results/A12B.sam ./results/A12B.bam ./results/A12B.sorte..
## V12A ./results/V12A.sam ./results/V12A.bam ./results/V12A.sorte..
## V12B ./results/V12B.sam ./results/V12B.bam ./results/V12B.sorte..
## samtools_index
## <character>
## M1A ./results/M1A.sorted..
## M1B ./results/M1B.sorted..
## A1A ./results/A1A.sorted..
## A1B ./results/A1B.sorted..
## V1A ./results/V1A.sorted..
## ... ...
## M12B ./results/M12B.sorte..
## A12A ./results/A12A.sorte..
## A12B ./results/A12B.sorte..
## V12A ./results/V12A.sorte..
## V12B ./results/V12B.sorte..In an ‘R-centric’ rather than a ‘CWL-centric’ workflow design the connectivity
among workflow steps is established by creating the downstream targets files
automatically (see Figure 3). Each step can uses the previous step outfiles as
an input, establishing connectivity among the steps in the workflow. By chaining
several SYSargsList steps together one can construct complex workflows involving
many sample-level input/output file operations with any combination of command-line or R-based
software. Also, systemPipeR provides features to automatically and systematically build this
connection, providing security that all the samples will be managed efficiently
and reproducibly.
Alternatively, a CWL-centric workflow design can be used that defines all/most workflow steps with CWL workflow and parameter files. Due to time and space restrictions, the CWL-centric approach is not covered by this tutorial.
The following step generate data frame containing important read alignment statistics such as the total number of reads in the FASTQ files, the number of total alignments, as well as the number of primary alignments in the corresponding BAM files.
This constructor function LineWise requires the step_name and the R code-based
under the code argument.
The R code should be enclosed by braces ({}) and separated by a new line.
appendStep(sal) <- LineWise(code = {
fqpaths <- getColumn(sal, step = "hisat2_mapping", "targetsWF", column = "FileName1")
bampaths <- getColumn(sal, step = "hisat2_mapping", "outfiles", column = "samtools_sort_bam")
read_statsDF <- alignStats(args = bampaths, fqpaths = fqpaths, pairEnd = TRUE)
write.table(read_statsDF, "results/alignStats.xls", row.names = FALSE, quote = FALSE,
sep = "\t")
}, step_name = "align_stats", dependency = "hisat2_mapping")Also, for printing and double-check the R code in the step, the codeLine method can be used:
codeLine(sal, "align_stats")## align_stats
## fqpaths <- getColumn(sal, step = "hisat2_mapping", "targetsWF", column = "FileName1")
## bampaths <- getColumn(sal, step = "hisat2_mapping", "outfiles", column = "samtools_sort_bam")
## read_statsDF <- alignStats(args = bampaths, fqpaths = fqpaths, pairEnd = TRUE)
## write.table(read_statsDF, "results/alignStats.xls", row.names = FALSE, quote = FALSE, sep = "\t")One interesting feature showed here is the getColumn method that allows
extracting the information for a workflow instance. Those files can be used in
an R code, as demonstrated below
getColumn(sal, step = "hisat2_mapping", "outfiles", column = "samtools_sort_bam")## M1A M1B
## "./results/M1A.sorted.bam" "./results/M1B.sorted.bam"
## A1A A1B
## "./results/A1A.sorted.bam" "./results/A1B.sorted.bam"
## V1A V1B
## "./results/V1A.sorted.bam" "./results/V1B.sorted.bam"
## M6A M6B
## "./results/M6A.sorted.bam" "./results/M6B.sorted.bam"
## A6A A6B
## "./results/A6A.sorted.bam" "./results/A6B.sorted.bam"
## V6A V6B
## "./results/V6A.sorted.bam" "./results/V6B.sorted.bam"
## M12A M12B
## "./results/M12A.sorted.bam" "./results/M12B.sorted.bam"
## A12A A12B
## "./results/A12A.sorted.bam" "./results/A12B.sorted.bam"
## V12A V12B
## "./results/V12A.sorted.bam" "./results/V12B.sorted.bam"0.3.2 Build workflow from R Markdown
The workflow can be created by importing the steps from an R Markdown file.
As demonstrated above, it is required to initialize the project with SPRproject function.
importWF function will scan and import all the R chunk from the R Markdown file
and build all the workflow instances. Then, each R chuck in the file will be
converted in a workflow step.
sal <- SPRproject()## Creating directory '/home/dcassol/z_tests/GEN242_toPush/content/en/tutorials/systempiper/rnaseq/.SPRproject'
## Creating file '/home/dcassol/z_tests/GEN242_toPush/content/en/tutorials/systempiper/rnaseq/.SPRproject/SYSargsList.yml'sal <- importWF(sal, file_path = "systemPipeRNAseq.Rmd")## Reading Rmd file##
## ---- Actions ----## Ignore none-R chunks at line: 25## Checking chunk eval values## Checking chunk SPR option## Ignore non-SPR chunks: 34, 43, 61, 96, 120, 187, 204, 297, 656, 684, 702, 710, 721## Parse chunk code## Now importing step 'load_SPR'
## Now importing step 'preprocessing'
## Now importing step 'trimming'
## Now importing step 'fastq_report'
## Now importing step 'hisat2_index'
## Now importing step 'hisat2_mapping'
## Now importing step 'align_stats'
## Now importing step 'bam_IGV'
## Now importing step 'create_db'
## Now importing step 'read_counting'
## Now importing step 'sample_tree'
## Now importing step 'run_edger'
## Now importing step 'custom_annot'
## Now importing step 'filter_degs'
## Now importing step 'venn_diagram'
## Now importing step 'get_go_annot'
## Now importing step 'go_enrich'
## Now importing step 'go_plot'
## Now importing step 'heatmap'
## Now importing step 'sessionInfo'Let’s explore the workflow to check the steps:
stepsWF(sal)
dependency(sal)
codeLine(sal)
targetsWF(sal)0.3.2.1 Third-party software tools
Current, systemPipeR provides the param file templates for third-party software tools.
A list is provided in the following table.
| Tool Name | Description | Step |
|---|---|---|
| bwa | BWA is a software package for mapping low-divergent sequences against a large reference genome, such as the human genome. | Alignment |
| Bowtie2 | Bowtie 2 is an ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences. | Alignment |
| FASTX-Toolkit | FASTX-Toolkit is a collection of command line tools for Short-Reads FASTA/FASTQ files preprocessing. | Read Preprocessing |
| TransRate | Transrate is software for de-novo transcriptome assembly quality analysis. | Quality |
| Gsnap | GSNAP is a genomic short-read nucleotide alignment program. | Alignment |
| Samtools | Samtools is a suite of programs for interacting with high-throughput sequencing data. | Post-processing |
| Trimmomatic | Trimmomatic is a flexible read trimming tool for Illumina NGS data. | Read Preprocessing |
| Rsubread | Rsubread is a Bioconductor software package that provides high-performance alignment and read counting functions for RNA-seq reads. | Alignment |
| Picard | Picard is a set of command line tools for manipulating high-throughput sequencing (HTS) data and formats such as SAM/BAM/CRAM and VCF. | Manipulating HTS data |
| Busco | BUSCO assesses genome assembly and annotation completeness with Benchmarking Universal Single-Copy Orthologs. | Quality |
| Hisat2 | HISAT2 is a fast and sensitive alignment program for mapping NGS reads (both DNA and RNA) to reference genomes. | Alignment |
| Tophat2 | TopHat is a fast splice junction mapper for RNA-Seq reads. | Alignment |
| GATK | Variant Discovery in High-Throughput Sequencing Data. | Variant Discovery |
| STAR | STAR is an ultrafast universal RNA-seq aligner. | Alignment |
| Trim_galore | Trim Galore is a wrapper around Cutadapt and FastQC to consistently apply adapter and quality trimming to FastQ files. | Read Preprocessing |
| TransDecoder | TransDecoder identifies candidate coding regions within transcript sequences. | Find Coding Regions |
| Trinity | Trinity assembles transcript sequences from Illumina RNA-Seq data. | denovo Transcriptome Assembly |
| Trinotate | Trinotate is a comprehensive annotation suite designed for automatic functional annotation of transcriptomes. | Transcriptome Functional Annotation |
| MACS2 | MACS2 identifies transcription factor binding sites in ChIP-seq data. | Peak calling |
| Kallisto | kallisto is a program for quantifying abundances of transcripts from RNA-Seq data. | Read counting |
| BCFtools | BCFtools is a program for variant calling and manipulating files in the Variant Call Format (VCF) and its binary counterpart BCF. | Variant Discovery |
| Bismark | Bismark is a program to map bisulfite treated sequencing reads to a genome of interest and perform methylation calls in a single step. | Bisulfite mapping |
| Fastqc | FastQC is a quality control tool for high throughput sequence data. | Quality |
| Blast | BLAST finds regions of similarity between biological sequences. | Blast |
Remember, if you desire to run any of these tools, make sure to have the respective software installed on your system and configure in the PATH. You can check as follows:
tryCL(command = "grep")0.4 How to run a Workflow
0.4.1 Setup and requirements
To go through this tutorial, you need the following software installed:
- R (version >=4.1.2)
- systemPipeR package (version >=2.0.8)
- Hisat2 (version >= 2.1.0)
If you desire to build your pipeline with any different software, make sure to have the respective software installed and available in your PATH. To make sure if the configuration is correct, on test it with:
tryCL(command = "hisat2") ## 'All set up, proceed!'0.4.2 Running the workflow
For running the workflow, runWF function will execute all the command lines
store in the workflow container.
sal <- runWF(sal)This essential function allows the user to choose one or multiple steps to be
executed using the steps argument. However, it is necessary to follow the
workflow dependency graph. If a selected step depends on a previous step(s) that
was not executed, the execution will fail.
sal <- runWF(sal, steps = c(1, 3))Also, it allows forcing the execution of the steps, even if the status of the
step is 'Success' and all the expected outfiles exists.
Another feature of the runWF function is ignoring all the warnings
and errors and running the workflow by the arguments warning.stop and
error.stop, respectively.
sal <- runWF(sal, force = TRUE, warning.stop = FALSE, error.stop = FALSE)0.4.2.1 Parallelization on clusters
This section of the tutorial provides an introduction to the usage of the
systemPipeR features on a cluster.
Alternatively, the computation can be greatly accelerated by processing many files in parallel using several compute nodes of a cluster, where a scheduling/queuing system is used for load balancing.
The resources list object provides the number of independent parallel cluster
processes defined under the Njobs element in the list. The following example
will run 18 processes in parallel using each 4 CPU cores.
If the resources available on a cluster allow running all 18 processes at the
same time, then the shown sample submission will utilize in a total of 72 CPU cores.
Note, runWF can be used with most queueing systems as it is based on utilities
from the batchtools package, which supports the use of template files (*.tmpl)
for defining the run parameters of different schedulers. To run the following
code, one needs to have both a conffile (see .batchtools.conf.R samples here)
and a template file (see *.tmpl samples here)
for the queueing available on a system. The following example uses the sample
conffile and template files for the Slurm scheduler provided by this package.
The resources can be appended when the step is generated, or it is possible to
add these resources later, as the following example using the addResources
function:
resources <- list(conffile=".batchtools.conf.R",
template="batchtools.slurm.tmpl",
Njobs=18,
walltime=120, ## minutes
ntasks=1,
ncpus=4,
memory=1024, ## Mb
partition = "short"
)
sal <- addResources(sal, c("hisat2_mapping"), resources = resources)
sal <- runWF(sal)Note: The example is submitting the jog to short partition. If you desire to
use a different partition, please adjust accordingly.
0.4.3 Visualize workflow
systemPipeR workflows instances can be visualized with the plotWF function.
This function will make a plot of selected workflow instance and the following information is displayed on the plot:
- Workflow structure (dependency graphs between different steps);
- Workflow step status, *e.g.* `Success`, `Error`, `Pending`, `Warnings`;
- Sample status and statistics;
- Workflow timing: running duration time. If no argument is provided, the basic plot will automatically detect width, height, layout, plot method, branches, etc.
plotWF(sal)0.4.4 Technical report
systemPipeR compiles all the workflow execution logs in one central location,
making it easier to check any standard output (stdout) or standard error
(stderr) for any command-line tools used on the workflow or the R code stdout.
Also, the workflow plot is appended at the beginning of the report, making it
easier to click on the respective step.
sal <- renderLogs(sal)0.4.5 Scientific report
systemPipeR compiles all the workflow execution logs in one central location,
making it easier to check any standard output (stdout) or standard error
(stderr) for any command-line tools used on the workflow or the R code stdout.
Also, the workflow plot is appended at the beginning of the report, making it
easier to click on the respective step.
sal <- renderLogs(sal)0.5 Workflow initialization with templates
Workflow templates are provided via systemPipeRdata and GitHub. Instances of these
workflows can be created with a single command.
0.5.1 RNA-Seq sample
Load the RNA-Seq sample workflow into your current working directory.
library(systemPipeRdata)
genWorkenvir(workflow = "rnaseq")
setwd("rnaseq")0.5.1.1 Create the workflow
This template provides some common steps for a RNAseq workflow. One can add, remove, modify
workflow steps by operating on the sal object.
sal <- SPRproject()
sal <- importWF(sal, file_path = "systemPipeRNAseq.Rmd", verbose = FALSE)Workflow includes following steps:
- Read preprocessing
- Quality filtering (trimming)
- FASTQ quality report
- Alignments:
HISAT2(or any other RNA-Seq aligner) - Alignment stats
- Read counting
- Sample-wise correlation analysis
- Analysis of differentially expressed genes (DEGs)
- GO term enrichment analysis
- Gene-wise clustering
0.5.1.2 Run workflow
sal <- runWF(sal)Workflow visualization
plotWF(sal)Report generation
sal <- renderReport(sal)
sal <- renderLogs(sal)0.5.2 ChIP-Seq sample
Load the ChIP-Seq sample workflow into your current working directory.
library(systemPipeRdata)
genWorkenvir(workflow = "chipseq")
setwd("chipseq")Workflow includes following steps:
- Read preprocessing
- Quality filtering (trimming)
- FASTQ quality report
- Alignments:
Bowtie2orrsubread - Alignment stats
- Peak calling:
MACS2 - Peak annotation with genomic context
- Differential binding analysis
- GO term enrichment analysis
- Motif analysis
0.5.2.1 Create the workflow
This template provides some common steps for a ChIPseq workflow. One can add, remove, modify
workflow steps by operating on the sal object.
sal <- SPRproject()
sal <- importWF(sal, file_path = "systemPipeChIPseq.Rmd", verbose = FALSE)0.5.2.2 Run workflow
sal <- runWF(sal)Workflow visualization
plotWF(sal)Report generation
sal <- renderReport(sal)
sal <- renderLogs(sal)0.5.3 VAR-Seq sample
Load the VAR-Seq sample workflow into your current working directory.
library(systemPipeRdata)
genWorkenvir(workflow = "varseq")
setwd("varseq")Workflow includes following steps:
- Read preprocessing
- Quality filtering (trimming)
- FASTQ quality report
- Alignments:
gsnap,bwa - Variant calling:
VariantTools,GATK,BCFtools - Variant filtering:
VariantToolsandVariantAnnotation - Variant annotation:
VariantAnnotation - Combine results from many samples
- Summary statistics of samples
0.5.3.1 Create the workflow
This template provides some common steps for a VARseq workflow. One can add, remove, modify
workflow steps by operating on the sal object.
sal <- SPRproject()
sal <- importWF(sal, file_path = "systemPipeVARseq.Rmd", verbose = FALSE)0.5.3.2 Run workflow
sal <- runWF(sal)Workflow visualization
plotWF(sal)Report generation
sal <- renderReport(sal)
sal <- renderLogs(sal)0.5.4 Ribo-Seq sample
Load the Ribo-Seq sample workflow into your current working directory.
library(systemPipeRdata)
genWorkenvir(workflow = "riboseq")
setwd("riboseq")Workflow includes following steps:
- Read preprocessing
- Adaptor trimming and quality filtering
- FASTQ quality report
- Alignments:
HISAT2(or any other RNA-Seq aligner) - Alignment stats
- Compute read distribution across genomic features
- Adding custom features to workflow (e.g. uORFs)
- Genomic read coverage along transcripts
- Read counting
- Sample-wise correlation analysis
- Analysis of differentially expressed genes (DEGs)
- GO term enrichment analysis
- Gene-wise clustering
- Differential ribosome binding (translational efficiency)
0.5.4.1 Create the workflow
This template provides some common steps for a RIBOseq workflow. One can add, remove, modify
workflow steps by operating on the sal object.
sal <- SPRproject()
sal <- importWF(sal, file_path = "systemPipeRIBOseq.Rmd", verbose = FALSE)0.5.4.2 Run workflow
sal <- runWF(sal)Workflow visualization
plotWF(sal, rstudio = TRUE)Report generation
sal <- renderReport(sal)
sal <- renderLogs(sal)0.6 Version information
Note: the most recent version of this tutorial can be found here.
sessionInfo()## R Under development (unstable) (2021-10-25 r81105)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /home/dcassol/src/R-devel/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] magrittr_2.0.3 batchtools_0.9.15
## [3] ape_5.6-2 ggplot2_3.3.5
## [5] systemPipeR_2.1.27 ShortRead_1.53.1
## [7] GenomicAlignments_1.31.2 SummarizedExperiment_1.25.3
## [9] Biobase_2.55.2 MatrixGenerics_1.7.0
## [11] matrixStats_0.62.0 BiocParallel_1.29.21
## [13] Rsamtools_2.11.0 Biostrings_2.63.3
## [15] XVector_0.35.0 GenomicRanges_1.47.6
## [17] GenomeInfoDb_1.31.7 IRanges_2.29.1
## [19] S4Vectors_0.33.17 BiocGenerics_0.41.2
## [21] BiocStyle_2.23.1
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-157 bitops_1.0-7 webshot_0.5.3
## [4] httr_1.4.2 RColorBrewer_1.1-3 progress_1.2.2
## [7] tools_4.2.0 backports_1.4.1 bslib_0.3.1
## [10] DT_0.22 utf8_1.2.2 R6_2.5.1
## [13] DBI_1.1.2 colorspace_2.0-3 withr_2.5.0
## [16] prettyunits_1.1.1 tidyselect_1.1.2 compiler_4.2.0
## [19] rvest_1.0.2 cli_3.2.0 formatR_1.12
## [22] xml2_1.3.3 DelayedArray_0.21.2 bookdown_0.26
## [25] sass_0.4.1 scales_1.2.0 checkmate_2.1.0
## [28] rappdirs_0.3.3 systemPipeRdata_1.99.10 systemfonts_1.0.4
## [31] stringr_1.4.0 digest_0.6.29 svglite_2.1.0
## [34] rmarkdown_2.13 jpeg_0.1-9 pkgconfig_2.0.3
## [37] htmltools_0.5.2 fastmap_1.1.0 htmlwidgets_1.5.4
## [40] rlang_1.0.2 rstudioapi_0.13 jquerylib_0.1.4
## [43] generics_0.1.2 hwriter_1.3.2.1 jsonlite_1.8.0
## [46] crosstalk_1.2.0 dplyr_1.0.8 RCurl_1.98-1.6
## [49] kableExtra_1.3.4 GenomeInfoDbData_1.2.8 Matrix_1.4-1
## [52] Rcpp_1.0.8.3 munsell_0.5.0 fansi_1.0.3
## [55] lifecycle_1.0.1 stringi_1.7.6 yaml_2.3.5
## [58] zlibbioc_1.41.0 grid_4.2.0 parallel_4.2.0
## [61] crayon_1.5.1 lattice_0.20-45 hms_1.1.1
## [64] knitr_1.38 pillar_1.7.0 base64url_1.4
## [67] codetools_0.2-18 glue_1.6.2 evaluate_0.15
## [70] latticeExtra_0.6-29 data.table_1.14.2 remotes_2.4.2
## [73] BiocManager_1.30.16 png_0.1-7 vctrs_0.4.1
## [76] gtable_0.3.0 purrr_0.3.4 assertthat_0.2.1
## [79] xfun_0.30 viridisLite_0.4.0 tibble_3.1.6
## [82] ellipsis_0.3.2 brew_1.0-70.7 Funding
This project is funded by NSF award ABI-1661152.